#############################################
#Robustness check: include voters in the same household as volunteers
#############################################

rm(list = ls())

#Functions
source("functions.R")

#Packages
installPackageNotFound("stargazer")
installPackageNotFound("lmtest")
installPackageNotFound("sandwich")

#############################################
#Read in pair data
#############################################
#Pairs
pair_df <- read.csv("voter_pairs.csv", stringsAsFactors = F)

#############################################
#Paired treatment effect
#############################################
#Model variables
vars <- c("BinaryTurnout", "Black", "Hispanic", "Asian", "Male", "Under30",
          "VotePropensity", "Partisanship")

#Data groups
data_list <- list("in_district" = which(pair_df$Volunteer_In_District == 1),
                  "out_district" = which(pair_df$Volunteer_In_District == 0))

#Specifications
itt_base <- paste0("BinaryTurnout_diff ~1")
itt_controls <- paste0("BinaryTurnout_diff ~",
                       paste0(vars[2:length(vars)], "_diff", collapse = "+"))
itt_assignment <- paste0("BinaryTurnout_diff~ -1 + factor(assignment_factor)")
itt_assignment_controls <- paste0(itt_assignment, "+",
                          paste0(vars[2:length(vars)], "_diff", collapse = "+"))
#List of formulas
formula_list <- c(itt_base,
                  itt_controls,
                  itt_assignment,
                  itt_assignment_controls)

#Empty list to hold model results
model_list <- list()
se_list <- list()

for(k in 1:length(data_list)){
  for(i in 1:length(formula_list)){
    #Subset to volunteer group
    this_group <- pair_df[data_list[[k]], ]

    #Calculate within-pair differences
    for(this_var in vars){
      this_group[, paste0(this_var, "_diff")] <-
        this_group[, paste0(this_var, "_treatment")] -
        this_group[, paste0(this_var, "_control")]
    }

    #Treatment effect at the pair level
    this_model <- lm(formula(formula_list[[i]]), data = this_group)

    #Robust SEs
    this_se <- coeftest(this_model, vcov = vcovHC(this_model, type = "HC3"))

    #Store results
    model_list[[names(data_list)[k]]][[i]] <- this_model
    se_list[[names(data_list)[k]]][[i]] <- this_se
  }
}

#############################################
#ITT Latex table
#############################################
title <- "Intention-to-treat Models"

#Full sample model
latex_output <- stargazer(model_list[["in_district"]][[1]],
model_list[["in_district"]][[2]],
model_list[["out_district"]][[1]],
model_list[["out_district"]][[2]],
se = list(se_list[["in_district"]][[1]][, "Std. Error"],
          se_list[["in_district"]][[2]][, "Std. Error"],
          se_list[["out_district"]][[1]][, "Std. Error"],
          se_list[["out_district"]][[2]][, "Std. Error"]
          ),
omit = paste0(vars[2:length(vars)], "_diff"),
covariate.labels = c("Intercept"),
omit.stat = c("ser", "rsq", "f"),
add.lines = list(c("Voter Covariates",
                   "\\multicolumn{1}{c}{N}",
                   "\\multicolumn{1}{c}{Y}",
                   "\\multicolumn{1}{c}{N}",
                   "\\multicolumn{1}{c}{Y}")),
column.separate = c(1,1,1,1),
dep.var.labels = "",
dep.var.caption = "",
notes.append = FALSE,
column.sep.width = "0pt",
no.space = TRUE,
title = title,
star.cutoffs = c(0.05, 0.01, 0.001),
type = "text")

#############################################
#CATE Latex table
#############################################
title <- "CATE Models"

#Full sample model
latex_output <- stargazer(model_list[["in_district"]][[3]],
                          model_list[["in_district"]][[4]],
                          model_list[["out_district"]][[3]],
                          model_list[["out_district"]][[4]],
                          se = list(se_list[["in_district"]][[3]][, "Std. Error"],
                                    se_list[["in_district"]][[4]][, "Std. Error"],
                                    se_list[["out_district"]][[3]][, "Std. Error"],
                                    se_list[["out_district"]][[4]][, "Std. Error"]
                          ),
                          omit = paste0(vars[2:length(vars)], "_diff"),
                          omit.stat = c("ser", "rsq", "f"),
                          covariate.labels = c("Order 1", "Order 2", "Order 3",
                                               "Order 4+"),
                          add.lines = list(c("Voter Covariates",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}")),
                          column.separate = c(1,1,1,1),
                          dep.var.labels = "",
                          dep.var.caption = "",
                          notes.append = FALSE,
                          column.sep.width = "0pt",
                          no.space = TRUE,
                          title = title,
                          star.cutoffs = c(0.05, 0.01, 0.001),
                          type = "text")
